Instalación de paquetes
Primero cargamos los paquetes usethis y devtools para descargar los paquetes desde github.
library(usethis)
library(devtools)Luego, usando la funcion install_github() pasando como arentos la dirección del repositorio en el cual está el paquete, y force = TRUE para evitar problemas durante la instalación.
install_github('https://github.com/sergmonic/AtomicWeights/',force = TRUE)
install_github('https://github.com/sergmonic/uLPAQ/',force = TRUE)Una vez instalados los paquetes, podemos agregarlos al entorno de trabajo utilizando la función library().
library(AtomicWeights)
library(uLPAQ)Podemos acceder a las páginas de ayuda para cada una de las funciones mediante la función help() o ?
help(package="AtomicWeights")
help(package="uLPAQ")Procedimiento para la evaluación y estimación según la GUM
Identificar el mesurando
Identificar las fuentes de incertidumbre y la posible existencia de fuentes dominantes.
Clasificar las fuentes de incertidumbre en tipo A o tipo B.
Expresar la relación matemática entre el mesurando \(Y\) y las variables de entrada \(X_i\) de las cuales depende:
\[\begin{align} Y = f(X_1, X_2, ..., x_n) \end{align}\]
Esta expresión debe contener cada cantidad que puede contribuir en la incertidumbre del resultado de la medición.
Estimar las incertidumbres asociadas a cada funte de error (contribuciones). Convertir a incertidumbres estándar.
Combinar las incertidumbres para obtener la incertidumbre combinada.
Para sumas:
\[\begin{align} y = a + b, y = a - b \end{align}\] \[\begin{align} u_y = \sqrt{u^2_a + u^2_b} \end{align}\] Para Multiplicaciones:
\[\begin{align} y = ab, y = a/b \end{align}\] \[\begin{align} u_y/y = \sqrt{(u_a/a)^2 + (u_b/b)^2} \end{align}\]
Calcular La Incertidumbre expandida (U) utilizando un factor de cubrimiento (k)
\[\begin{align} U = ku \end{align}\]
Modelo Matematico del mesurando e identificación de fuentes de incertidumbre.
Mesurando: Concentración molar de NaOH: obtenida mediante titulación volumétrica empleando KHP y con seguimiento potenciométrico.
Fuentes:
- Masa de KHP empleada para la estandarización (
m_kHP) en g - Masa molar de KHP (
W_KHP) en g/mol - Volumen gastado de solución de NaOH (
V_NaOH) en mL - Pureza del KHP (
P_KHP) de forma fraccional (g/g) - Repetibilidad (R) obtenida de X muestras independientes (fraccional/relativa)
- Masa de KHP empleada para la estandarización (
Para relacionar las fuentes de incertidumbre con la incertidumbre del mesurando, podemos dibujar un diagrama Ishikawa de causa y efecto, o espina de pescado.
Para construir este diagrama en R, podemos utilizar la función cause.and.effect del paquete qcc:
library(qcc)
cause.and.effect(cause=list(m_KHP = '(g)',
W_KHP = '(g/mol)',
V_NaOH = c('material volumetrico',
'punto final'),
P_KHP = '(g/g)',
Repetibilidad= ''),
effect = 'Incertidumbre NaOH')Modelo matemático
Para el cálculo de la concentración \(C_{NaOH}\) de NaOH de una solución, la cual se valora haciendo reaccionar \(m_{KHP}\) g de biftalato de potasio, el cual tiene una pureza de \(m_{KHP}\) con un volumen de \(V_{NaOH}\) de la solución problema:
\[\begin{align} C_{NaOH} = \frac{1000 \; \cdot \; m_{KHP} \; \cdot \; P_{KHP} \; \cdot \ R}{W_{KHP} \; \cdot \; V_{NaOH} } \end{align}\]
Podemos crear la función que relaciona la concentración del hidróxido de sodio C_NaOH con el factor de conversión, utilizando la función function() creando como argumentos cada uno de los parámetros del factor de conversión:
C_NaOH<-function(m_KHP,
W_KHP,
V_NaOH,
P_KHP,
R){
1000*m_KHP*P_KHP*R/(W_KHP*V_NaOH)
}Incertidumbre en masas molares
Fórmula molecular de KHP: \(KC_8O_4H_5\)
Para determinar la masa molar de un compuesto, y la incertidumbre asociada a esta masa molar, podemos utilizar las funciones getAtomicWeight y getMolarMass, disponibles en el paquete AtomicWeights.
Para revisar más información sobre dichas funciones, podemos revisar las siguientes páginas de ayuda:
help("getAtomicWeight")
help("getMolarMass")W_KHP=getMolarMass(atomsType = c("K", "C", "O", "H"),
c(1,8,4,5))Volumen de NaOH gastado en la valoración
Incertidumbre en mediciones de volumen
Para las mediciones de volumen en material volumétrico de vidrio, podemos utilizar la función instVolumetrico. Para acceder a la página de ayuda de esta función podemos ejecutar el siguiente comando:
help("instVolumetrico")Incertidumbre la determinación del punto final de la titulación (volumen de equivalencia)
Para determinar la incertidumbre en la determinación del punto final de la titulación, e.g. en una valoración potenciométrica con electrodo de pH, el punto de inflexión en la curva pH vs volumen de agente titulante, podemos utilizar la función puntoFinal del paquete uLPAQ.
Esta función determina el punto final mediante tres algoritmos distintos, y devuelve una lista con tres tablas, o data.frame. Cada una correspondiente a un algoritmo de determinación del punto de inflexión:
d1numcontiene el resultado en \(x\) (volumen de titulante) determinado como el máximo de la primera derivada numérica hacia adelanterLogcontiene el resultado del punto de inflexión empleando el modelo de regresión logística disponible en el paquetenplredecontiene el punto final de la titulacion estimado mediante el estimador de distancia máxima, disponible en el paqueteinflection
PF=puntoFinal(datosAnalisis$valoracionNaOH$V,
datosAnalisis$valoracionNaOH$pH,
intervalo_x = c(10,13))## value
## npar 5
## params.bottom 2.49021e-02
## params.top 9.33177e-01
## params.xmid 1.01057e+01
## params.scal 2.28388e+00
## params.s 1.72950e+01
## GOF 9.90708e-01
## weightedGOF 9.99698e-01
## stdErr 0.03920632
## weighted stdErr 0.02549547
## trapezoid 2.13022
## Simpson 2.11183
## xInfl 10.6478
## yInfl 0.368472
## Log10(IC50) 1.030466
## IC50 1.07267e+01
## [95%] [1.06867e+01 | 1.07697e+01]
## date (Y-m-d) 2021-06-28
## nplr version 0.1.7 (2016-12-25)
## R version 4.0.5 (2021-03-31)
PF## $d1num
## PE u_PE
## 1 10.65 0.04082483
##
## $rLog
## PE u_PE
## 1 10.6478 0.01846104
##
## $ede
## PE u_PE
## 1 10.7 0.3819993
Incertidumbre por el uso de material volumétrico (bureta de capacidad de 25 mL).
para determinar la incertidumbre por medir el volumen determinado con la función de puntoFinal en una bureta, podemos utilizar la función instVolumetrico, como se muestra a continuación:
bureta_NaOH <- instVolumetrico(nominal=25,
subdivision =0.05 ,
tipo="bureta",
clase="A_AS",
volumen=PF$d1num$PE,
temp = 17.5)obtención de V_NaOH considerando material volumétrico y obtención de PF
Una vez tenemos el aporte a la incertidumbre por las dos fuentes para el volumen de titulante gastado en la valoración:
Incertidumbre asociada a la determinación del punto final, contenida en
PF$d1num$u_PEIncertidumbre asociada a el uso de material volumetrico contenida en
bureta_NaOH$u_V.
Ahora podemos calcular la incertidumbre asociada al volumen como incertidumbre combinada de las dos fuentes principales, punto final y medición en bureta.
\[\begin{align} u_{V_{NaOH}} = \sqrt{u_{PF}^2 + u_{bureta}^2} \end{align}\]
.V_NaOH=PF$d1num$PE
.u_V_NaOH=sqrt(PF$d1num$u_PE^2+bureta_NaOH$u_V^2)
V_NaOH=data.frame(V_NaOH=.V_NaOH,
u_V_NaOH=.u_V_NaOH)
V_NaOH## V_NaOH u_V_NaOH
## 1 10.65 0.04274444
Ejemplos con otros tipos de material volumétrico
Si se llega a dar el caso en el que en el procedimiento experimental se utilizan otros tipos de material de vidrio para realizar mediciones de volumen, podemos utilizar la función instVolumetrico. A continuación algunos ejemplos.
balon100 <- instVolumetrico(nominal=100,
tipo="balon",
clase="A",
temp=16)pipeta10 <- instVolumetrico(nominal=10,
tipo="pipeta",
clase="A_AS",
temp=16)micropipeta <- instVolumetrico(nominal=1,
tipo="micropipeta",
volumen=1,
temp=16)Incertidumbre en medición de masa de KHP: mediciones gravimétricas
Dentro del paquete uLPAQ tenemos disponibles las funciones u_masa y masa_minima relacionadas con el uso de balanzas analíticas.
masa_minimanos sirve para determinar cual es la mínima cantidad de masa que podemos pesar en una balanza con determinadas propiedades para cumplir con la condición de un valor máximo de incertidumbre relativa.u_masanos sirve para determinar la incertidumbre asociada a medir la masa de un sólido que tiene una densidad asociada, y usando una balanza analítica con propiedades determinadas.
Las balanzas que se pueden utilizar en este paquete, se encuentran en la tabla balanzasLPAQ dentro del paquete uLPAQ:
| Característica | ME204 | AL54 | AES220 | XPE204 | XPE205 | Unidad |
|---|---|---|---|---|---|---|
| Capacidad máxima | 2.2e+02 | 5.1e+01 | 2.2e+02 | 2.2e+02 | 2.2e+02 | g |
| Legibilidad | 1.0e-01 | 1.0e-01 | 1.0e-01 | 1.0e-01 | 1.0e-02 | mg |
| Repetibilidad | 1.0e-01 | 1.0e-01 | 1.0e-01 | 7.0e-02 | 3.0e-02 | mg |
| Desviacio´n de la linealidad | 2.0e-01 | 2.0e-01 | 3.0e-01 | 2.0e-01 | 1.0e-01 | mg |
| Sensibilidad de la deriva te´rmica | 2.0e-06 | 2.5e-06 | 2.0e-06 | 1.0e-06 | 1.0e-06 | g/(g °C) |
| Tolerancia en la Sensibilidad | 1.5e-06 | 1.5e-06 | 1.5e-06 | 1.0e-06 | 1.0e-06 | g/g |
En donde:
\(Capacidad\;Máxima\) es la masa máxima que se puede pesar en la balanza, reportada por el fabricante
\(Legibilidad\) es la mínima división de escala, \(d\), reportada por el fabricante.
\(Repetibilidad\) Asociada con la dispersión de mediciones repetidas de masas patrón, expresada por el fabricante o determinada en una calibración.
La \(desviación \; de\;la \;linealidad \;\) nos da información sobre que tanto varía la relación lineal entre la masa nominal y la lectura de la balanza.
\(la \;deriva \;en \;la \;sensibilidad \;por\; la\; deriva\; termica\) nos habla de como cambia la pendiente de la relación masa nominal \(\rightarrow\) lectura, por un cambio en la temperatura.
la \(Tolerancia\; en \;la \;Sensibilidad\) expresa un intervalo dentro del cual pueden estar los valores de la sensibilidad, con una probabilidad dada.
A partir de estas propiedades, es que las funciones masa_minima y u_masa pueden calcular masas en funcion de incertidumbres relativas, e incertidumbres y masas corregidas.
Para saber más información sobre masa_minima:
help(masa_minima)Por ejemplo, podemos calcular la cantidad de masa mínima de biftalato de potasio, con una densidad de \(1636 \frac{Kg}{m^3}\), para cumplir con un máximo de incertidumbre relativa para la medición de masa \(0.1\;g\):
masa_minima(u_r = 0.1,densidad = 1636,balanza = "ME204")A continuación, calculamos la incertidumbre asociada a la medición de 0.1026\;g en la balanza ME204, de un sólido con densidad de \(1636 \frac{Kg}{m^3}\) (biftalato de potasio):
m_KHP=u_masa(lectura = 0.1026,
densidad = 1636,
balanza = "ME204")La función u_masa nos devuelve una lista de tablas, teniendo en la primera posición de esta lista, un data.frame con dos columnas, en la primera columna la masa corregida por el efecto de flotación, y en la segunda columna la incertidumbre asociada a la medición de esta masa, en esta balanza.
str(m_KHP)## List of 3
## $ m_s:'data.frame': 1 obs. of 2 variables:
## ..$ ms : num 0.103
## ..$ u_ms: num 0.000255
## $ Bu :'data.frame': 2 obs. of 5 variables:
## ..$ Bu : num [1:2] 1.00 2.09e-05
## ..$ d_s : num [1:2] 1636 47.2
## ..$ d_as: num [1:2] 0.8847 0.0224
## ..$ d_ac: num [1:2] 1.2 0.0129
## ..$ d_pc: num [1:2] 8006 10
## $ Ws :'data.frame': 1 obs. of 7 variables:
## ..$ ws : num 103
## ..$ u_ws : num 0.255
## ..$ u_Res: num 0.0408
## ..$ u_rep: num 0.1
## ..$ u_NL : num 0.231
## ..$ u_TS : num 8.89e-11
## ..$ u_ETS: num 5.33e-10
Incertidumbre en la pureza del biftalato
La función pureza devuelve estimados de puereza en incertidumbre basados en especificaciones técnicas del reactivo. Debemos especificar los valores como fracciones de masa, en el intervalo \(0 \leq pureza \leq 1\).
help("pureza")P_KHP = pureza(0.999)
P_KHP## P u_P
## 1 0.9997071 0.0002357023
Incertidumbre por repetibilidad
Para realizar el cálculo de la concentración de la solución de hidróxido de sodio, sin ser afectado por el valor de la repetibilidad, definimos la repetibilidad como R =1.
Dado que en el experimento de este ejercicio solo se realizó una determinación, sin réplicas, la incertidumbre asociada a la repetiblidad se define como u_R=0.
En el caso en el que se tengan réplicas, podemos reportar la incertidumbre asociada a la repetiblidad tal como se reporta en la literatura:
\[\begin{align} u_R = \frac{s}{\sqrt{n}} \end{align}\]
en donde \(s\) es la desviación estándar de las réplicas, (el resultado de concentración) y \(n\) es el número de réplicas.
R=1
u_R=0.0Estimación de la incertidumbre combinada
Una vez hemos determinado todas las fuentes de incertidumbre y estimado las incertidumbres asociadas a dichas mediciones podemos estimarlas de la siquiente manera:
- de forma general
\[\begin{align} u^2_c (y) = \sum_{i=1}^N \left(\frac{\delta f}{\delta x_i}\right)^2 \cdot u^2(x_i) \end{align}\]
\[\begin{align} u_c (y) = \left( \sum_{i=1}^N \left(\frac{\delta f}{\delta x_i}\right)^2 \cdot u^2(x_i) \right)^\frac{1}{2} \end{align}\]
\[\begin{align} u_c (y) = \sqrt{ \sum_{i=1}^N \left(\frac{\delta f}{\delta x_i}\right)^2 \cdot u^2(x_i) } \end{align}\]
En donde \(\left(\frac{\delta f}{\delta x_i}\right)\) es el coeficiente de sensiblidad de \(f\) respecto a \(x_i\).
- En este caso
\[\begin{align} y = f(x_1, x_2, x_3,...,x_n) = C_{NaOH} (m_{KHP}, P_{KHP}, R, W_{KHP}, V_{NaOH}) \end{align}\]
\[\begin{align} u^2_c (C_{NaOH}) = \sum_{i=1}^5 \left(\frac{\delta f}{\delta x_i}\right)^2 \cdot u^2(x_i) \end{align}\]
Creación de una función que estima los coeficientes de sensibilidad para una función
Para encontrar como cambia la concentración de hidróxido de sodio, como respuesta a los cambios en cada una de las variables de infuencia: m_KHP,W_KHP,V_NaOH,P_KHPyR Podemos utilizar la función Deriv del paquete con el mismo nombre, pasando como primer argumento la función que contiene el modelo matemático en R, C_NaOH, y como segundo argumento un vector de caracteres, que contiene los nombres de las variables con respecto a las cuales la función C_NaOH debe ser diferenciada.
library(Deriv)
c_i=Deriv(C_NaOH,
c("m_KHP","W_KHP","V_NaOH","P_KHP","R"))al imprimir el contenido de la funcion c_i vemos una descripcion de los pasos que se toman para calcular las derivadas de la función con una información de entrada determinada (los valores que toman las variables m_KHP,W_KHP,V_NaOH,P_KHPyR)
c_i ## function (m_KHP, W_KHP, V_NaOH, P_KHP, R)
## {
## .e1 <- V_NaOH * W_KHP
## .e2 <- m_KHP * P_KHP
## .e3 <- .e1^2
## .e4 <- .e2 * R
## c(m_KHP = 1000 * (P_KHP * R/.e1), W_KHP = -(1000 * (.e4 *
## V_NaOH/.e3)), V_NaOH = -(1000 * (.e4 * W_KHP/.e3)), P_KHP = 1000 *
## (m_KHP * R/.e1), R = 1000 * (.e2/.e1))
## }
Valores de los coeficientes de sensibilidad, obtenido a partir de los valores de la variables de influencia
Al tener lista la función c_i, podemos darle los valores de las variables de influencia para este experimento, dado que ya las definimos dentro de R, e incluso corregimos algunas como la pureza y la masa del biftalato en pasos anteriores. En este caso c_i a partir de los valores definidos calculará la magnitud de cada uno de los coeficientes de sensibilidad y lo guardará en cada posición del vector v_c_i:
v_c_i <- c_i(m_KHP$m_s$ms,W_KHP$MM,V_NaOH$V_NaOH,P_KHP$P,R)
v_c_i## m_KHP W_KHP V_NaOH P_KHP R
## 0.4596462032 -0.0002310156 -0.0044298726 0.0471919652 0.0471781430
vector recogiendo las incertdumbres estándar de la diferentes fuentes
Ya tenemos calculadas las incertidumbres de cada fuente, pero éstas están guardadas en cada objeto que creamos en los pasos anteriores, para juntarlas, podemos llamarlas por su nombre y guardarlas en un vector u_i en el mismo orden en el que están en el vector v_c_i.
u_i=c(m_KHP$m_s$u_ms,W_KHP$u_MM,V_NaOH$u_V_NaOH,P_KHP$u_P,u_R)
u_i## [1] 0.0002550597 0.0047143761 0.0427444424 0.0002357023 0.0000000000
Incertidumbre combinada
Una vez tenemos los valores de los coeficientes de sensiblidad, y de las incertidumbres, podemos combinarlas utilizando la ecuación:
\[\begin{equation} u_c (C_{NaOH}) = \sqrt{\left( \frac{\delta C_{NaOH}}{\delta m_{KHP}}\right)^2 \cdot u^2(m_{KHP}) + \left( \frac{\delta C_{NaOH}}{\delta P_{KHP}}\right)^2 \cdot u^2( P_{KHP}) + \left( \frac{\delta C_{NaOH}}{\delta R}\right)^2 \cdot u^2( R) + \\ \left( \frac{\delta C_{NaOH}}{\delta W_{KHP}}\right)^2 \cdot u^2(W_{KHP}) + \left( \frac{\delta C_{NaOH}}{\delta V_{NaOH}}\right)^2 \cdot u^2(V_{NaOH}) } \end{equation}\]
u_c=sqrt(sum((v_c_i*u_i)^2))
u_c## [1] 0.0002229884
Obtención de incertidumbre expandida
Una medida adicional de incertidumbre que cumple con el requerimiento de proporcionar un intervalo con nuestro resultado (para reportar la concentración) es la incertidumbre expandida , denotada como \(U\). Ésta se obtiene al multiplicar la incertidumbre combinada, \(u_c(C_{NaOH})\) por un factor de \(k\):
\[\begin{align} U = k \cdot u_c(C_{NaOH}) \end{align}\]
k=2 # Se expande asumiendo una distribución normal, para una probabilidad de confianza del 95 %
U=k*u_c
U # Incertidumbre expandida## [1] 0.0004459768
Contribuciones de cada fuente
Una vez tenemos la incertidumbre combinada y el valor de cada una de las fuentes de incertidumbre, podemos calcular la fracción relativa de cada una de las fuentes:
\[\begin{align} contribución \;(x_i) = \frac{\left( \frac{\delta C_{NaOH}}{\delta x_i}\right)^2 \cdot u(x_i)^2}{u_c(C_{NaOH})} \end{align}\]
contribuciones=(v_c_i*u_i)^2/u_c^2
pie(contribuciones,main="Contribuciones de fuentes de incertidumbre")barplot(contribuciones*100,
xlab="Fuentes de incertidumbre",
ylab="Aporte a incertudumbre combinada (%)",
main="Contribuciones de fuentes de incertidumbre")Valor del mensurando
Una vez hemos determinado la incertidumbre expandida, podemos reportar el valor de concentración utilizando la función c_NaOH que contiene el modelo matemático que relaciona la concentración de la solución, con las variables de influencia m_KHP,W_KHP,V_NaOH,P_KHPyR.
.C_NaOH=C_NaOH(m_KHP$m_s$ms,W_KHP$MM,V_NaOH$V_NaOH,P_KHP$P,R)mensurando_C_NaOH=round(data.frame(C_NaOH=.C_NaOH,U_C_NaOH=U),5)
mensurando_C_NaOH## C_NaOH U_C_NaOH
## 1 0.04718 0.00045
Ejemplo con datos reales
Para este ejemplo se utilizaran los datos resultantes de valoraciones repetidas de una muestra de cascara de huevo. Para determinar el porcentaje de carbonato de calcio \((CaCO_3)\) se determinó el contenido de calcio por gravimetría y por volumetría. En total, 20 grupos de laboratorio determinaron el contenido de calcio, mediante distintas metodologias de la siguiente manera:
- los grupos del 1 al 10 tomaron 100 mg de la muestra en la balanza
ME204 - los grupos del 11 al 20 grupos tomaron 200 mg de la muestra en la balanza
AES220
Cada muestra se llevó a volumen en balones tipo A de 50mL.
Luego de esto, se tomaron alicuotas de la siguiente manera:
- Todos los grupos tomaron alicuota de 10 mL en una pipeta tipo A.
Luego se llevó a cabo la valoración de cada una de las alícuotas de la siguiente manera:
Las muestras 1 a 10 se llevaron a pH 10 y se titularon con EDTA 0.100 \(\pm\) 0.001 M (k=2) en una bureta de 25 \(\pm\) 0.05 mL tipo A
Las muestras 11 a 20 precipitaron con oxalato de amonio, filtraron, secaron, calcinaron a 1000 ºC y se pesaron en la balanza
ME204.
Un resumen de las variables de procedimiento y los resultados experimentales se encuentran en la siguiente tabla:
| Masa muestra (g) | balanza inicial | Alícuota (mL) | Clase pipeta alícuota | Capacidad de la bureta (mL) | balanza final | Volumen EDTA (mL) | Masa final (g) | replica |
|---|---|---|---|---|---|---|---|---|
| 0.5 | ME204 | 10 | A | 25 | NA | 8.75 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.59 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.99 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.28 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.67 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.70 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.16 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.18 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.17 | NA | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.01 | NA | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2617 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.1612 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2614 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2623 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2632 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2621 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2617 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2613 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2614 | 1 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2646 | 1 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.77 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.42 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.82 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.25 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.47 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.38 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.15 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.93 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.07 | NA | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.77 | NA | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2631 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.1601 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2605 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2620 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2622 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2611 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2606 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2600 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2609 | 2 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2641 | 2 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.92 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.42 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.78 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.10 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.33 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.42 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.73 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.88 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 9.03 | NA | 3 |
| 0.5 | ME204 | 10 | A | 25 | NA | 8.84 | NA | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2599 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.1606 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2603 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2613 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2630 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2627 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2633 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2605 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2630 | 3 |
| 2.5 | AES220 | 10 | A | NA | ME204 | NA | 0.2634 | 3 |
Podemos importar esta tabla de datos, disponible en este enlace:
library(downloader)
url1 <- 'https://github.com/pipoelpipas/Statistics-in-R-2/raw/main/CLT-incertidumbre-clase/tablaEjemplo.csv'
download(url1, destfile = 'tablaEjemplo.csv')
datos <- read.csv('tablaEjemplo.csv', sep = ',', header = T)Una vez tenemos cargados los datos, podemos dividirlos en grupos mas pequeños, por ejemplo un grupo de volumetría y otro de gravimetría:
library(dplyr)
gravimetria <- filter(datos, balanza.final == 'ME204')
volumetria <- filter(datos, Capacidad.de.la.bureta..mL. == 25)Analisis exploratorio
Podemos realizar análisis exploratorio de datos, tal como se ha revisado en el curso a cada uno de los grupos de experimentos, antes de proceder a calcular el contenido de carbonato de calcio en la cáscara de huevo.
masas <- unlist(select(gravimetria, Masa.final..g.))
volumenes <- unlist(select(volumetria, Volumen.EDTA..mL.))
par(mfrow = c(2,2))
boxplot(masas, main = 'boxplot de masas')
hist(masas)
boxplot(volumenes, main = 'boxplot de volumenes')
hist(volumenes)En este análisis gráfico exploratorio, podemos ver que la distribución de las masas para las titulaciones gravimétricas de aleja de la normalidad debido a datos atípicos, podemos preguntarle a R cuales son, primero calculando el valor crítico inferior:
summary(masas)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1601 0.2605 0.2616 0.2518 0.2629 0.2646
cuartiles <- quantile(masas)
cuartiles## 0% 25% 50% 75% 100%
## 0.160100 0.260525 0.261550 0.262925 0.264600
IQR <- cuartiles[4]-cuartiles[2]
critico.inferior<- cuartiles[2]-(1.5*IQR)
which(masas < critico.inferior)## Masa.final..g.2 Masa.final..g.12 Masa.final..g.22
## 2 12 22
Ahora sabemos que las masas del grupo 12 (el segundo de gravimetría) son datos atípicos, así que podemos rechazarlos. para extraer de nuestro set de datos estos resultados podemos usar el siguiente código:
gravimetria <- gravimetria[-c(2,12,22),]Ahora podemos volver a evaluar las distribuciones empíricas de los experimentos:
masas <- unlist(select(gravimetria, Masa.final..g.))
par(mfrow = c(2,2))
boxplot(masas, main = 'boxplot de masas')
hist(masas)
boxplot(volumenes, main = 'boxplot de volumenes')
hist(volumenes) Ahora podemos ver que no hay datos atípicos sesgando la distribución de los resultados experimentales.
Una vez realizamos este análisis exploratorio, podemos proceder a calcular los promedios de las réplicas y sus respectivas desviaciones estándar (Éstas serán útiles para la incertidumbre dada por repetibilidad para cada experimento)
- Para las volumetrias:
volumenes1 <- filter(volumetria, replica ==1) %>% select(Volumen.EDTA..mL.)%>% unlist()
volumenes2 <- filter(volumetria, replica ==2) %>% select(Volumen.EDTA..mL.)%>% unlist()
volumenes3 <- filter(volumetria, replica ==3) %>% select(Volumen.EDTA..mL.)%>% unlist()
volumenes.promedio <- (volumenes1+volumenes2+volumenes3)/3
sd.volumen <- vector('numeric',10)
for(i in 1:10){
sd.volumen[i] <- sd(c(volumenes1[i],volumenes2[i],volumenes3[i]))
}
volumenes.promedio <- as.data.frame(volumenes.promedio)
# sd.volumen
# volumenes.promedio$sd.replicas <- sd.volumen- Para las gravimetrías:
masas1 <- filter(gravimetria, replica ==1) %>% select(Masa.final..g.)%>% unlist()
masas2 <- filter(gravimetria, replica ==2) %>% select(Masa.final..g.)%>% unlist()
masas3 <- filter(gravimetria, replica ==3) %>% select(Masa.final..g.)%>% unlist()
masas.promedio <- (masas1+masas2+masas3)/3
sd.masa <- vector('numeric',10)
for(i in 1:10){
sd.masa[i] <- sd(c(masas1[i],masas2[i],masas3[i]))
}
sd.masa## [1] 0.0016041613 0.0005859465 0.0005131601 0.0005291503 0.0008082904
## [6] 0.0013576941 0.0006557439 0.0010969655 0.0006027714 NA
Una vez tenemos separados nuestros datos, podemos proceder a estimar la incertidumbre teniendo encuenta los pasos recomendados por la GUM
Volumetría
Modelo Matematico del mesurando e identificación de fuentes de incertidumbre
Para el caso de la volumetria, podemos expresar el modelo de la siguiente forma:
\[\begin{align} \%CaCO_3 = \frac{0.1 \;\cdot\; FD_{muestra} \;\cdot\; V_{EDTA}\;\cdot\; W_{CaCO_3}\;\cdot\;C_{EDTA} \;\cdot\;R}{m_{muestra}} \end{align}\]
Podemos entonces crear esta función en R:
CaCO3_p.p <- function(FD,
V_EDTA,
W_CaCO3,
C_EDTA,
R,
m_muestra){
(0.1*FD*V_EDTA*W_CaCO3*C_EDTA*R)/(m_muestra)
} En donde:
\(FD\_{muestra}:\) es el factor de dilución correspondiente a la toma de la alícuota con una pipeta de 10 mL a partir de un volumen de aforo de 50 mL en un balón volumétrico.
\(W\_CaCO_3:\) Es la masa molar del carbonato de calcio
\(V\_EDTA:\) es el volumen de EDTA gastado en la titulación, utilizando un material volumétrico (bureta) y un método para determinar el punto final (indicador, o mediante una curva potenciométrica)
\(C_{EDTA}:\) La concentración de la solución de EDTA valorada con una incertidumbre expandida de 0.001 M (k = 2)
\(m\_muestra:\) es la masa inicial de la muestra, medida en una balanza específica.
\(Repetibilidad:\) La repetibilidad de los experimentos expresada como la desviación de las réplicas, con un valor de 1 y un aporte a la incertidumbre de \(U_R = \dfrac{s}{\sqrt{n}}\).
Incertidumbre por masa molar de \(CaCO_3\)
Para determinar la incertidumbre asociada a la masa molar de carbonato de calcio, podemos utilizar de nuevo la función getMolarMass del paquete AtomicWeights.
library(AtomicWeights)
W_CaCO3 <- getMolarMass(atomsType = c('Ca','C','O'),
c(1,1,3)
)
W_CaCO3## MM u_MM
## 1 100.0868 0.004091947
Volumen de EDTA gastado en la valoración.
Incertidumbre en la medición del punto final
En este caso, tenemos dos variaciones, en el laboratorio solo había un electrodo de ion selectivo para calcio, entonces solo el primer grupo pudo hacer la determinación mediante titulación potenciométrica por triplicado.
El resto de grupos, realizaron la determinación del punto final mediante indicador, como se muestra en la siguiente imagen:
Entonces, podemos:
Determinar los puntos finales con su respectiva incertidumbre para las replicas del primer grupo.
Dados los puntos finales de los resultados de potenciometría, calcular la incertidumbre dada por la bureta.
Una vez tengamos estos valores, podemos hallar las incertidumbres de los puntos determinados con indicador, tal como se ha reportado en literatura.
Para ilustrar de nuevo el caso de la determinación de la incertidumbre en el punto final de la titulación, podemos cargar los datos contenidos en el archivo curvasEDTA.csv disponible en este link
Primero cargamos los datos a R:
library(downloader)
url2 <- 'https://github.com/pipoelpipas/Statistics-in-R-2/raw/main/CLT-incertidumbre-clase/curvasEDTA.csv'
download(url2, destfile = 'curvasEDTA.csv')
curvasEDTA <- read.csv('curvasEDTA.csv', sep = ',', header = TRUE)Podemos graficar las curvas de titulación:
plot(curvasEDTA[,1],
curvasEDTA[,2],
type = 'l',
xlab = 'Volumen de EDTA (mL)',
ylab = 'pCa',
xlim = c(5,10),
col = 1,
ylim = c(1,9)
)
par(new =T)
plot(curvasEDTA[,3],
curvasEDTA[,4],
type = 'l',
xlab = '',
ylab = '',
xlim = c(5,10),
col = 2,
ylim = c(1,9)
)
par(new =T)
plot(curvasEDTA[,5],
curvasEDTA[,6],
type = 'l',
xlab = '',
ylab = '',
xlim = c(5,10),
col = 3,
ylim = c(1,9)
)Luego, podemos utilizar la función puntoFinal del paquete uLPAQ para estimar la incertidumbre por la determinación del punto final a partir de esta curva:
library(uLPAQ)
puntos.finales <- list(3)
par(mfrow=c(1,3))
for(i in c(1,3,5)){
puntos.finales[i] <- puntoFinal(curvasEDTA[,i],curvasEDTA[,i+1])
}puntos.finalesUna vez tenemos el volumen dado por la potenciometría, podemos promediar los resultados, y estimar la incertidumbre de este promedio:
volumenes.promedio$volumenes.promedio[1] <- mean(puntos.finales[[1]]$PE,puntos.finales[[3]]$PE,puntos.finales[[5]]$PE)
volumenes.promedio <- as.data.frame(volumenes.promedio)
volumenes.promedio$U.PF.Class[1] <- sqrt(puntos.finales[[1]]$u_PE+puntos.finales[[3]]$u_PE+puntos.finales[[5]]$u_PE)Para determinar las incertidumbres de las determinaciones mediante uso del indicador, se pueden utilizar dos procedimientos:
- El procedimiento clásico, el cual se utiliza cuando No hay un valor de referencia mediante el cual se determina la incertidumbre en la determinación del punto final mediante:
\[\begin{align} u_{Qpf} = \frac{R}{2\sqrt{3}} \end{align}\]
En donde \(U_{Qpf}\) es la incertidumbre en el punto final y \(R\) corresponde a la resolución del istrumento volumétrico. En la literatura se encuentra reportado que este procedimiento no aplica para todos los casos y solo aplica si el error sistematico es igual a la resolución del instrumento.
- El procedimiento mediante el cual se hace uso de un valor de referencia
\[\begin{align} u_{Qpf} = \frac{D}{\sqrt{3}} \end{align}\]
Donde D corresponde a la estimación del error sistemático debido al sistema de detección, y no al instrumento de medición del agente titulante.
Entonces, para los últimos 9 grupos, podemos determinar la incertidumbre de las dos maneras, ya que tenemos un valor de referencia con las curvas de potenciometría.
volumenes.promedio$U.PF.Class[2:10] <- 0.05/2*sqrt(3)A su vez podemos calcular las incertidumbres con el método de referencia:
volumenes.promedio$U.PF.Ref[1] <- volumenes.promedio$U.PF.Class[1]
volumenes.promedio$U.PF.Ref[2:10] <- abs(volumenes.promedio$volumenes.promedio[1]-volumenes.promedio$volumenes.promedio[2:10])/sqrt(3)
par(mfrow=c(1,2))
barplot(volumenes.promedio$U.PF.Class, ylim = c(0,0.4))
barplot(volumenes.promedio$U.PF.Ref, ylim = c(0,0.4))Incertidumbre por el uso de material volumétrico (bureta de capacidad de 25mL)
library(uLPAQ)
U_buretas <- vector('list',10)
for(i in 1:10){
U_buretas[[i]] <- instVolumetrico(nominal=25,
subdivision =0.05 ,
tipo="bureta",
clase="A_AS",
volumen=volumenes.promedio$volumenes.promedio[i],
temp = 17.5)$u_V
}
U_buretas <- unlist(U_buretas)
for(i in 1:10){
volumenes.promedio$u_bureta[i] <- as.numeric(U_buretas[i])
}Incertidumbre en el volumen de agente titulante: EDTA teniendo en cuenta la determinación del punto final y el uso de material volumétrico:
Una vez tenemos las incertidumbres para todos los grupos, podemos calcular la incertidumbre combinada para la determinación del volumen de titulante:
u_V_EDTA <- vector('numeric',10)
u_V_EDTA <- sqrt(volumenes.promedio$U.PF.Ref^2+volumenes.promedio$u_bureta^2)
volumenes.promedio$u_V_EDTA <- u_V_EDTA
par(mfrow = c(1,3))
barplot(volumenes.promedio$U.PF.Ref, main = 'u punto final' , ylim = c(0,0.5))
barplot(volumenes.promedio$u_bureta, main = 'u bureta', ylim = c(0,0.5))
barplot(volumenes.promedio$u_V_EDTA, main = 'u volumen EDTA' , ylim = c(0,0.5))Incertidumbre en mediciones de masa: masa inicial de muestra de cáscara de huevo.
Ya que todos los grupos tuvieron la destreza de pesar 0.500 g, podemos calcular la incertidumbre dada por la pesada de esta muestra de la siguiente manera:
- Para este experimento se determino la denisdad de la cascara de huevo como: $ 1085 $
masas.corregidas <- vector('numeric',10)
u_masas <- vector('numeric',10)
for(i in 1:10){
masas.corregidas[i] <- u_masa(lectura=0.500,
densidad = 1085,
balanza = "ME204")$m_s$ms
}
for(i in 1:10){
u_masas[i] <- u_masa(lectura=0.500,
densidad = 1085,
balanza = "ME204")$m_s$u_ms
}
volumenes.promedio$masas_corregidas <- masas.corregidas
volumenes.promedio$u_masas <- u_masasIncertidumbre por dilución y toma de alícuota
Para este procedimiento se utilizaron:
- un balón de 50 \(\pm\) 0.050 mL clase A
- una pipeta de 10 \(\pm\) 0.02 mL clase A
balon50 <- instVolumetrico(nominal =50,
tipo="balon",
clase="A",
temp=17.5)
pipeta10 <- instVolumetrico(nominal=10,tipo="pipeta",clase="A_AS",temp=17.5)
u_dilucion <- sqrt(balon50$u_V^2 + pipeta10$u_V^2)
volumenes.promedio$u_dilucion <- rep(u_dilucion,10)Incertidumbre por concentración de la solución de EDTA
La concentración de la solución estandarizada reporta una incertidumbre expandida de: \(\pm\) 0.001M con un k =2
podemos calcular la incertidumbre estandar:
u_C_EDTA <- 0.001/2
volumenes.promedio$u_C_EDTA <- rep(u_C_EDTA,10)Incertidumbre por repetibilidad
Debido a que cada experimento se realizo por triplicado, podemos usar la incertidumbre asociada a la repetibilidad para cada experimento como:
\(\frac{s}{\sqrt{n}}\)
Dado que ya estimamos las desviaciones estandar de las réplicas para cada grupo, podemos añadir a nuestro data.frame llamado volumenes.promedio en el cual estamos guardando por columnas, las incertidumbres de cada fuente para los 1 grupos que realizaron una volumetría.
volumenes.promedio$u_R <- sd.volumen/sqrt(10)Creación de la función que estima los coeficientes de sensibilidad
Dado el modelo matematico construido a partir del factor de conversión necesario para calcular el porcentaje de \(CaCO_3\), almacenado en la función CaCO3_p.p, podemos utilizar la función Deriv del paquete que lleva el mismo nombre para calcular las derivadas parciales del porcentaje de carbonato de calcio respecto a las variables de entrada:
library(Deriv)
c_i = Deriv(CaCO3_p.p,c('FD',
'V_EDTA',
'W_CaCO3',
'C_EDTA',
' R',
'm_muestra'))
c_i## function (FD, V_EDTA, W_CaCO3, C_EDTA, R, m_muestra)
## {
## .e2 <- C_EDTA * FD * R
## .e3 <- .e2 * V_EDTA
## c(FD = 0.1 * (C_EDTA * R * V_EDTA * W_CaCO3/m_muestra), V_EDTA = 0.1 *
## (.e2 * W_CaCO3/m_muestra), W_CaCO3 = 0.1 * (.e3/m_muestra),
## C_EDTA = 0.1 * (FD * R * V_EDTA * W_CaCO3/m_muestra),
## ` R` = 0, m_muestra = -(0.1 * (.e3 * W_CaCO3/m_muestra^2)))
## }
Valores de los coeficientes de sensibilidad, obtenido a partir de los valores de la variables de influencia
coef.matrix <- matrix(nrow = 10, ncol = 6)
for(i in 1:10){
coef.matrix[i,] <- c_i(FD = 50/10,
V_EDTA = volumenes.promedio$volumenes.promedio[i],
W_CaCO3 = W_CaCO3$MM,
C_EDTA = 0.100,
R = 1,
m_muestra = 0.500)
}Matriz recogiendo las incertidumbres estándar de las distintas fuentes
u_matrix = data.frame(FD = volumenes.promedio$u_dilucion,
V_EDTA = volumenes.promedio$u_V_EDTA,
W_CaCO3 = rep(W_CaCO3$u_MM,10),
C_EDTA = volumenes.promedio$u_C_EDTA,
R = volumenes.promedio$u_R,
m_muestra = volumenes.promedio$u_masas)Ahora podemos calcular un vector de incertidumbres combinadas para cada uno de los grupos de laboratorio que desarrollaron la valoración mediante volumetría.
Primero, creamos un vector para guardar las incertidumbres:
u_c <- vector('numeric',10)Luego calculamos fila por fila, la incertidumbre combinada:
for(i in 1:10){
u_c[i] <- sqrt(sum((coef.matrix[i,]*u_matrix[i,])^2))
}Incertidumbres expandidas
Una vez tenemos el valor de las incertidumbres combinadas para cada muestra, podemos escoger un factor de cobertura de 2 y expresar la incertidumbre expandida.
k=2
U=k*u_cContribuciones de cada fuente
contribuciones.vol <- matrix(nrow = 10, ncol = 6)
for(i in 1:10){
for(j in 1:6){
contribuciones.vol[i,j] <- ((coef.matrix[i,j]*u_matrix[i,j])^2)/u_c[i]^2
}
}
colnames(contribuciones.vol) <- colnames(u_matrix)Con el calculo de todas las fuentes de incertidumbre, y la incertidumbre combinada, podemos graficar y compararlas entre muestras:
#win.graph()
#par(mfrow = c(5,2))
par(las=2)
par(mfrow = c(2,5))
for(i in 1:10){
barplot(contribuciones.vol[i,]*100,
xlab="(%)",
ylab=" ",
main= paste('contribuciones muestra #', as.character(i)),
horiz =T)
}Reporte de las concentraciones de cada grupo:
primero, calculamos las concentraciones de cada uno de las muestras, a partir de los volumenes promediados:
concentraciones <- vector('numeric', 10)
for (i in 1:10){
concentraciones[i] <- CaCO3_p.p(FD = 50/10,
V_EDTA = volumenes.promedio$volumenes.promedio[i],
W_CaCO3 = W_CaCO3$MM,
C_EDTA =0.100,
R = 1,
m_muestra = 0.5
)
}mesurando.volumen <- data.frame(`% p/p CaCO3` = concentraciones,
`U` = U)
colnames(mesurando.volumen) <- c('CaCO3 (% p/p)','U')knitr::kable(mesurando.volumen, align = 'c')| CaCO3 (% p/p) | U |
|---|---|
| 87.57683 | 1.434347 |
| 93.51443 | 7.015757 |
| 94.58203 | 8.227625 |
| 90.04476 | 3.190305 |
| 93.14745 | 6.600974 |
| 93.04736 | 6.488053 |
| 91.97977 | 5.290721 |
| 91.61278 | 4.883224 |
| 92.64701 | 6.037378 |
| 91.57942 | 4.846317 |
Gravimetria
Construcción del modelo matemático e identificacion de fuentes de incertidumbre
En este caso el modelo matemático para el porcentaje de \(CaCO_3\) es:
\[\begin{align} \% \; CaCO_3 = \frac{m_f \;\cdot\; W_{CaCO_3}\;\cdot\;FD\;\cdot\;100}{W_{CaO}\;\cdot\; m_{muestra}} \end{align}\]
Podemos entonces en este caso realizar el análisis mediante el diagrama de causa y efecto de Ishikawa: